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Parameterization of subgrid-scale 
stress by the velocity gradient tensor 


By T, S. Lund AND E. A. Novikov 1 


1. Motivation and objectives 

In large eddy simulation (LES), the large scale motions are directly computed 
and the effects of the small scales are modeled. The term to be modeled is the 
subgrid-scale stress tensor, 

Tij = TZTUJ- UjU;, (1) 

which arises when the Navier-Stokes equations are spatially filtered (denoted by 
overbar) to remove the small scale information. Most large eddy simulations make 
use of the Smagorinsky (1963) eddy-viscosity model, 

T ij ~ \r kk 6 tj = -2(CA 2 |5|)5 0 , (2) 

where C is a non-dimensional constant, A is the grid spacing, and Sij is the resolved 
strain-rate. 

Although the Smagorinsky model has been in use for nearly thirty years, for 
roughly half that period it has been known that the model provides only a crude 
estimate for the stresses. This fact was first demonstrated by Clark et al. (1979), 
where direct numerical simulation (DNS) data for homogeneous isotropic turbulence 
was used to evaluate model predictions. Clark et al. found a correlation coefficient 
of approximately 0.2 when comparing predictions of the Smagorinsky model with 
the exact stresses. McMillan et al (1979) found that the correlation coefficient was 
even lower in homogeneous shear flow, being close to 0.1. Later, Piomelli et al. 
(1988) found similar results in turbulent channel flow. 

When contemplating these extremely low correlations, it may seem striking that 
the Smagorinsky model is successful when used in a large eddy simulation. The rea- 
son for the seemingly unwarranted accuracy is that, by construction, the Smagorin- 
sky model insures a net drain of energy from the large scales to the subgrid-scale 
motions. This is the primary objective of the subgrid-scale model, and as long as 
this requirement is met, reasonable results are evidently obtained. On the other 
hand, the Smagorinsky model provides poor predictions of the individual elements 
of the stress tensor. It is natural to expect that superior results could be obtained 
with a model that predicts the stress tensor more accurately. The aim of this work 
is to seek out potentially more accurate models. 

The Smagorinsky model is based on a molecular transport analogy where the 
stress is proportional to the rate of strain. The molecular analogy is a rather crude 
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model for turbulent transport, however, and it may be expected that a more accu- 
rate model for the subgrid-scale stress could be obtained if additional information 
were included. An obvious quantity to consider is the vorticity (or equivalently the 
rotation rate tensor). While molecular transport is unaffected by rotation, there is 
no reason to exclude rotation from a model for turbulent transport. Indeed, vortex 
stretching is believed to be the dominant mechanism by which turbulence transfers 
energy from large to smaller scales (Tennekes & Lumley 1971). 

-"The objective of this work is to construct and evaluate subgrid-scale models that 
depend on both the strain rate and the vorticity. This will be accomplished by first 
assuming that the subgrid- scale stress is a function of the strain and rotation rate 
tensors. Extensions of the Caley-Hamilton theorem can then be used to write the 
assumed functional dependence explicitly in the form of a tensor polynomial involv- 
ing products of the strain and rotation rates. Finally, use of this explicit expression 
as a subgrid-scale model will be evaluated using direct numerical simulation data 
for homogeneous, isotropic turbulence. 

2* Accomplishments 

2.1 Subgrid-scale stress as a tensor function of strain and rotation rates 

It is assumed that the residual stress may be written as a function of the strain 
and rotation rate tensors, as well as the unit isotropic tensor, viz. 


T i j — 

(3) 

where the strain and rotation rate tensors are defined as 


Ifdiii duj\ 
* ij ~ 2 y dx j + dXi) ’ 

(4a) 

1 (dui dwA 

Kii - 2 \dx } dxj' 

(46) 

and where u, is the resolved velocity. Note that S and R 

are the symmetric and 


antisymmetric parts of the velocity gradient tensor. Thus Eq. (3) can be alterna- 
tively interpreted as a relationship between the subgrid-scale stress and the velocity 
gradient tensor. 

For simplicity in the following development, a matrix notation is introduced. Let 
r, S, R, and I be the matrices associated with the corresponding tensor quantities 
(I is the identity matrix associated with 6{j). Tensor contractions correspond to 
matrix multiplications in the following way: 

SR = SikRkj* R 2 = RikRkj , etc., (5a) 

tr(SR 2 ) — Sij RjkRki, etc. (5b) 

The most general expression for Eq. (3) is an infinite tensor polynominal con- 
taining terms of the form S ai R^S 0rj R^ 2 . . where and ft are positive integers. 
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Each term in the series is multiplied by a coefficient that may be a function of 
the simultaneous invariants of S and R (discussed below). The Caley-Hamilton 
theorem of matrix algebra (see Pipes & Hovanessian 1969) states that terms in the 
polynomial beyond a certain order axe redundant, and thus the infinite polynomial 
is equivalent to one involving a finite number of terms. Spencer and Rivlin (1959) 
have determined the surviving terms when one tensor is expressed in terms of two 
others, as done here. Since the subgrid-scale stress is symmetric, only the sym- 
metric tensors in this list need to be considered. When the Results of Spencer and 
Rivlin are applied to the present situation and the results properly symmetrized, 
the following tensors and associated invariants arise: 


mi = S, 

m 2 = S 2 , 

m 3 = R 2 , 

m 4 = SR — RS, 

m 5 = S 2 R — RS 2 , 

m 6 = I, 

m- =SR 2 + R 2 S, 

mg = RSR 2 — R 2 SR, 

mg = SRS 2 - S Z RS, 

mio = S 2 R 2 + R 2 S 2 , 

m„ = RS 2 R 2 -R z S 2 R, 



h = tr (S 2 ) , 
h = tr (S 3 ) , 
h = tr (S 2 R 2 ) , 


I 2 =tr(R 2 ), 
h = tr (SR 2 ) , 
h = tr (S 2 R 2 SR) . 


( 7 ) 


These results were presented earlier by Pope (1975) in connection with models for 
the Reynolds averaged Navier-Stokes equations. The invariant / 6 , however, was 
not included in Pope’s analysis. This invariant is anomalous in the sense that it is 
related to the other 5 but with an ambiguity in sign. The relationship is 


h = ±1 (4/j + 7 2 )7 3 - hl\ - I 3 + 47,(7 3 7 4 - 27,7 5 )7 2 + 
7i (/ 5 - hh)ll +57, 7 2 7 2 - 37 2 7 3 7 4 7 5 ) l '\ 


where I\ = /i/2, I 2 = /2/2, and 1 3 = 1$/ 3. The sign ambiguity can be resolved 
by replacing the ± above by the multiplicative factor sign(u?itl?2^3), where izq, W 2 , 
and w$ are the vorticity components expressed in the principal coordinate system 
of the strain rate tensor. (The function sign(x) returns either —1 for x < 0 or 1 
for x > 0.) Since resolution of the sign ambiguity requires information that is not 
contained in the invariants ij, I 2 , ... /s, the invariant Ie may be considered to be 
independent of the other 5. For this reason, we shall include 1$ in the subsequent 
analysis. 

The set of tensors displayed in Eq. (6) are complete in the sense that any symmet- 
ric polynomial involving products of S and R can be written as a linear combination 
of the 11 tensors, with the scalar multipliers expressed as polynomials of the 6 in- 
variants. The tensors are also independent in the sense that none of the 11 tensors 
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may be written as a linear combination of the other 10 if the scalar multipliers are 
restricted to be polynomials of the 6 invariants. If this restriction is relaxed slightly 
so that the scalar multipliers may be ratios of polynomials of the invariants, then 
under the conditions discussed below only 6 of the above 11 tensors are independent 
(see Rivlin and Ericksen (1955) for more details). To see this, consider expressing 
one of the 11 tensors as a linear combination of 6 others: 

m* = CiBiij * = 1,2,. ..6, k > 6, (8) 

where the tensors are ordered in any desired way, not necessarily as in Eq. (6). 
Due to symmetry, each of the tensors mi have only 6 unique elements. Thus Eq. 
(8) represents 6 algebraic equations for the 6 unknown coefficients C t . The solution 
can be written as 

Ci - [tr(niimj)] -1 tr(m*m>). (9) 

A unique solution will exist provided the above matrix of traces is non-singular, 
that is 

det[tr (miirij)) ^4 0. (10) 

Note that if Eq. (9) is solved by Cramer’s rule, the Ci will be expressed as a ratios 
of polynomials of the invariants listed in Eq. (7). 

Equation (10) can be violated under two conditions: when S has a repeated 
eigenvalue, or when two components of the vorticity, expressed in the principal 
coordinates of S, vanish. The first condition corresponds to an axisymmetric state 
of strain. The second corresponds to a situation where the rotation is confined to a 
single axis, and this axis is aligned with one of the principal directions of the strain 
rate. Although either of these conditions could be realized in a turbulent field, the 
probability of exactly satisfying either of them is rather remote. Indeed, when the 
tensor expansion was evaluated using direct numerical simulation data as described 
in the following section, the conditions for lack of independence were never satisfied 
exactly, even for 128 3 realizations. 

Assuming that Eq. (10) is satisfied, only the first 6 terms in Eq. (6) need to be 
considered. For incompressible flows, it is customary to model only the deviatoric 
part of t and combine the isotropic part with the pressure. We shall follow the 
precedent here and subtract the trace from each of the first 6 tensors in Eq. (6). 
The 6th term vanishes when this is done, leaving only the first 5 as a basis. This 
is consistent with the fact that a trace-free symmetric tensor has only 5 unique 
elements. This result could have been obtained equivalently by subtracting the 
trace from each of the 11 tensors at the outset and then showing that any tensor 
in the list can be written as a linear combination of the first 5. In any event, the 
stress can be written as 

t* =Ci A 2 |S|S + C 2 A 2 (S 2 )* + C 3 A 2 (R 2 )*+ 

C 4 A 2 (SR - RS) + C 5 A 2 ^|(S 2 R - RS 2 ), 


( 11 ) 
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where A is the grid spacing, |S| = y/ir (S 2 ), and ()* indicates the trace-free part. 
Use of the strain rate magnitude as a scaling factor was chosen somewhat arbitrarily. 
In theory, this is not an issue since the Ci can depend on all of the invariants in 
Eq. (7) and the correct scaling will be obtained if the Ci are written as functions 
of the invariants. In practice, it is difficult to find the dependence of the C t on the 
invariants, and thus the choice of the scaling becomes relevant. Several alternate 
scalings were tested and the results appeared to be quite insensitive to the particular 
choice. 

Since the expansion coefficients C, are non-dimensional, they can depend only on 
non-dimensional groupings of the invariants listed in Eq. (7). These are taken to 
be 


•si 


*2 


•S3 


S 4 


S 5 


tr(S 3 ) 
tr (S 2 ) 3/2 ' 
tr (R 2 ) 

tr(S 2 )’ 

tr (SR 2 ) 

tr (S 2 ) 1/2 tr(R 2 )' 
tr (S 2 R 2 ) 
tr (S 2 ) tr (R 2 ) 
tr (S 2 R 2 SR) 
[tr(S 2 )tr(R 2 )]<3/2)‘ 


(12a) 

(126) 

(12c) 

(12 d) 
(12c) 


2.2 Evaluation of the proposed model 

Equation (11) is an exact result that will hold as long a the basic assumption 
that the stress is expressible solely as a function of the strain and rotation rates is 
correct( i.e. Eq. (3)). Thus under this assumption, the stress r can be represented 
exactly in terms of the strain and rotation rate tensors, provided the coefficients 
Ci are known functions of the invariants listed in Eq. (12). The functional form of 
the dependence on the invariants is unknown, however, and can not be determined 
easily. If the assumption that the stress is expressible solely as a function of the 
strain and rotation rates is not correct, then the coefficients will be functions of the 
unknown quantities on which the stress really depends. In either case, it can be 
anticipated that it will be difficult to predict the spatial variation of the expansion 
coefficients. 

In the context of modeling, the expansion coefficients would most likely be as- 
signed constant values that reflect an overall “best-fit” for all points in the field. If 
the true coefficient values do not vary greatly in space, then taking them to be con- 
stant will be a reasonable approximation and a good representation of the stresses 
can be expected. On the other hand, if the coefficients vary greatly in space, then 
taking them to be constant would be a poor approximation and the model would be 
of little value. Thus, in practical terms, the utility of this approach depends on the 
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degree to which the expansion coefficients vary or, alternately, how well Eq. (11) is 
satisfied for fixed coefficients. 

The issue of coefficient variability was investigated through the use of direct 
numerical simulation (DNS) data of homogeneous, isotropic turbulence. By filtering 
the DNS field with a spectral cutoff filter, the subgrid-scale stress, as well as the 
“resolved” strain and rotation rates were computed exactly. The accuracy of Eq. 
(11) was then measured in two alternate ways: (1) by determining the expansion 
coefficients exactly at each grid point and then measuring their spatial variation and 
(2) by measuring the degree to which Eq. (11) was satisfied when the coefficients 
were assigned constant values. 

The homogeneous, isotropic data was generated with a pseudo-spectral code (Ro- 
gallo, (1981) on a 128 3 mesh. The energy spectrum was initialized according to 

m - h (I) exp H) ■ (13) 

This spectrum has its energy peak at wavenumber 8. The initial phases were chosen 
randomly, but in such a way that the divergence-free condition was satisfied (see 
Rogallo, (1981) for more details on the initial conditions). The flow was allowed 
to evolve freely for 2.9 small scale eddy turnover times, £ where A is the Taylor 
microscale and u* is the rms turbulence intensity, both based on the final field. Over 
this period of time, the total turbulent kinetic energy decayed by 33%. The final 
Taylor microscale Reynolds number (u'A/i/) was 45.3, and the velocity derivative 
skewness was -0.32. The final energy spectrum, scaled in Kolmogorov units, is 
plotted in Figure 1. 

Also shown in Figure 1 are the experimental data of Comte-Bellot and Corrsin 
(1971). The simulation results fit well with the experimental data. The tail-up in 
the simulated spectrum at high wavenumbers is a characteristic of spectral methods 
and is more pronounced when the dissipation range is not fully resolved, as in this 
case. It is generally felt (Rogallo, private communication) that the tail-up at high 
wavenumber will not adversely affect the data in the central portion of the spectrum 
used here. 

Following the procedure of Clark et a/.(!975), a synthetic LES velocity field was 
generated from the DNS data by filtering out the small scale motions. The filtering 
was achieved via spherical truncation in wave space where three quarters of the 
active high frequency modes were removed. The LES field thus corresponded to an 
isotropic simulation performed on a 128/ 4 — 32 cubed mesh. Denoting the filtering 
operation with an overbar, the subgrid-scale stress was determined by performing 
the operations in the de-aliased definition commonly used in spectral calculations, 

Tij - U~U~ - UiUy (14) 

The large scale strain and rotation rate tensors defined in Eq. (4) were determined 
by applying spectral derivative operators to the LES velocity field. 
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Wavenumber, kr) 

FIGURE 1. Energy spectrum from the DNS data base. The experimental data are 
taken from Comte-Bellot and Corrsin (1971). The vertical line corresponds to the 
scale at which the velocity field was filtered to generate the synthetic LES field. 


2.3 Analysis for variable coefficients 

With the subgrid-scale stress and the large scale strain and rotation rate tensors 
known, the expansion coefficients in Eq. (11) could be determined at each point in 
the field. This was done using a least-squares fitting procedure so that solutions 
could be obtained when less than all five of the tensors on the right hand side of 
Eq. (11) were used. To derive the least-squares expression, consider the error in 
satisfying Eq. (11) when an arbitrary number of tensors are used: 

E = Ciirii - r, (15) 

where m, are the model tensors on the right hand side of Eq. (11), and i = 1,2, ...n; 
1 < re < 5. The square of the error will be minimized with respect to the C x if the 
following condition is enforced 

~QQ. ix (E 2 ) = 0- (16) 

This condition leads to the following algebraic system for the coefficients: 

Ci = [tr(m,m J )) _ 1 tr(m > T). (17) 


If less than 5 model tensors are used, the subgrid-scale stress can not be represented 
exactly by Eq. (11). The following quantities are global measures of the error in 
this case: 
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Cr 



< tr ((t - M) 2 ) > 
< tr(r 2 ) > 


(19) 


where M = Cjirij is the composite model tensor and <> denotes a volume average. 
The quantity rj is the correlation coefficient between the exact and modeled stress, 
while e r is the rms error in the subgrid-scale stress. It is easy to show that rj and 
e r are related via e r = y'l - r? 2 . The variability of the coefficients themselves was 
measured in terms of the ratio of rms to mean value: 


C 


rms — 


V<c 2 > - <c > 2 
< c > 


( 20 ) 


2.3.1 Results 

As a first step, each of the 5 terms in Eq. (11) was considered separately. Mea- 
sures of the error as well as coefficient variability were recorded for each term. Next, 
the 10 possible pairs of terms in Eq. (11) were investigated. The 10 possible triplets, 
the 5 possible quadruplets, and finally all 5 terms together were tested. For each of 
the groupings, the combinations that resulted in the highest as well as the lowest 
correlation coefficients were selected for further study. Figure 2 shows these corre- 
lation coefficients as a function of the number of terms in the group. As expected, 
the correlation coefficient rises as more terms Eire added and is unity if all five terms 
are present. The differences between the best and worst correlation coefficients Eire 
rather slight, indicating that none of the terms are neither far superior nor far in- 
ferior to the rest. The terms forming the best and worst subsets are listed in Table 
1 . 


Number of terms 

Best combination 

Worst combination 

i 

i 

3 

2 

1,4 

2,3 

3 

1,4,5 

3, 4,5 

4 

1,2, 4,5 

2, 3, 4, 5 

5 

1, 2, 3, 4, 5 

1, 2, 3, 4, 5 


TABLE 1. Best and worst subsets of the model terms in Eq. (11) - variable 
coefficients. 

Although the differences between correlation coefficients obtained with the best 
and worst groupings are slight, there are some consistent trends if the terms are 
ranked by their relative importance. The best single term is term 1, which cor- 
responds to the Smagorinsky model. This term is present in each of the optimal 
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Number of model tensors 

FIGURE 2 . Correlation coefficients for the best and worst subsets of the terms in 
Eq. ( 11 ) - variable coefficients. 

groupings. The worst single term is the rotation rate squared (term 3 ). This term is 
also the last one to enter in the optimal groupings. Most of the intermediate terms 
follow the general trend that if they enter the optimal groupings when n terms are 
present, then they enter the worst grouping when 5 — n terms are used. 

Although ranking the various combinations of terms by their correlation coef- 
ficients is interesting, the more important issue is spatial variability of the corre- 
sponding expansion coefficients. The ratio of the coefficient rms to the coefficient 
mean for the optimal groupings of terms is shown in Figure 3 . It is clear that a 
substantial variation in each coefficient is required to achieve the least-squares fit. If 
only one term is included, the coefficient variation is roughly three times the mean. 
As more terms are included, the variation rapidly increases. When all five terms 
are included, the coefficient variation is enormous, ranging from roughly 10 times 
the mean for C\ to over 500 times the mean for C3. 

The rather large coefficient variation could in part be due to the the neglected 
dependence on the invariants listed in Eq. (12). It is conceivable that if this 
dependence were taken into account, the coefficient variability could be reduced. 
This issue is explored in the following section. 

2.4 Dependence on the invariants 

Each of the expansion coefficients in Eq. (11) can, in principal, depend of the five 
invariants listed in Eq. (12). Determining such a dependence in a five parameter 
space is a difficult task, however. One way to do this would be to divide the range 
of each invariant into m intervals, thereby partitioning the parameter space into m 5 
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Number of model tensors 


FIGURE 3. Coefficient variability for optimal subsets of the terms in Eq. (11). 

hypercubes. Using the DNS data, the coefficients could then be sorted according 
to their associated invariant values and the resulting sample averaged over each 
hypercube. Unfortunately, this procedure requires an enormous amount of data 
if reliable statistics are to be obtained. As an illustration, consider the following 
example. If 16 intervals are chosen for the discretization and a 128 3 DNS data base is 
used, there will be on the average only 128 3 /16 5 = 2 samples within each hypercube. 
This sample is clearly too small to provide meaningful statistics. Furthermore, it 
may be anticipated that many of the cubes will contain no data at all. If a larger 
data base or more realizations are used and if the number of intervals is reduced, it 
may be possible to obtain reliable statistics. If this were done, the difficult task of 
fitting the statistical data with some sort of multidimensional function would still 
remain. 

In view of these difficulties, two alternate approaches have been adopted here. 
In the first, the number of invariants was reduced to one by assuming that the 
stress depended only on the strain rate. The local smoothing procedure described 
above was used in this case since the sample size was large and the resulting one- 
dimensional function could be easily curve fit. In the second approach, the full 
problem was considered and a sophisticated regression algorithm was used to find 
any dependence as well as its associated functional form. 

2-4- 1 Dependence on a single invariant 

The question of dependence on the invariants can be answered completely if it 
is first assumed that the stress depends only on the strain rate. In this case, the 
Caley-Hamilton theorem states that the stress can be explicitly written as a linear 
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combination of S, S 2 , and I (all higher powers of S are related to these three terms). 
If only the deviatoric part of r is to be modeled, then only S and (S 2 )* are required. 
Thus Eq. (11) reduces to first two terms in this case. The corresponding coefficients, 
Ci and C 2 , could depend at most on the invariant sj listed in Eq. (12). 


0.50 


0.25 


-0.25 



2.50 

1.25 

0 

-1.25 



FIGURE 4. Scatter plot of model coefficients versus the invariant s\. Only every 
8 th data point in each direction is plotted. 

This abbreviated model was tested as in section 2.3, with the coefficients deter- 
mined locally. The resulting coefficients are plotted as a function of the associated 
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FIGURE 5. Model coefficient conditionally averaged on invariant S] . 


invariant value in Figure 4. There does not seem to be any strong dependence on 
the invariant in either case. If the raw data are averaged over narrow intervals in 
s, t however, a weak trend emerges. The results of such an averaging are shown in 
Figure 5. It appears that both Ci and C *2 depend linearly on the invariant sj. Also 
shown in Figure 5 is the rms fluctuation of the data about the smoothed curve as 
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well as the rms fluctuation about the coefficient mean value. There is no visible 
reduction in the fluctuation level if the dependence on the invariant is accounted 
for. Thus while a clear dependence on the invariant was found, it accounts for vary 
little of the coefficient variation. 

2.4*2 Dependence on all five invariants: Projection pursuit regression 

Although the local smoothing procedure described in Section 2.4 does not seem 
suitable for this problem, there are other numerical methods that can perform 
multi-variable regression, even for a moderate sample size. Perhaps the best of 
these methods is the projection pursuit regression algorithm developed by Friedman 
and Stuetzle in 1981. The algorithm consists of a numerical optimization routine 
that finds one dimensional projections of the original independent variables for 
which the best correlations with the dependent variable can be obtained. The 
dependent variable is then written as a sum of empirically determined functions 
of these projections. The method is quite robust and has been able to determine 
nonlinear relationships within a 5 parameter space using only 213 realizations (see 
Friedman and Stuetzle (1982) for more details and examples). The method has also 
been used by Meneveau et al. (1992) to search DNS data for improved subgrid-scale 
model parameterizations. 

The projection pursuit algorithm was used to search for coefficient dependence 
on the five invariants. The DNS data was used to determine the five expansion 
coefficients and five invariants at each point in the field. This data was then input 
to the projection pursuit algorithm and each coefficient was analyzed independently. 
For each coefficient, the numerical optimization routine was able to find projections 
for which the variance was minimized, but the reduction in variance never exceeded 
2%. Furthermore, the empirically determined functions of the invariants did not 
appear to have recognizable structure and contained many oscillations. This type of 
behavior is often an indication that the algorithm has only found a local minimum 
in field of noise. 

The results of the projection pursuit regression are consistent with the results 
presented in the previous section where the dependence on a single invariant was 
investigated. In both cases, the coefficients do appear to depend on the invariants, 
but that dependence is extremely weak. Accounting for this weak dependence on 
the invariants does not significantly reduce the coefficient spatial variation and thus 
is probably not worth pursuing further. More importantly, the large coefficient 
variation does not appear to be related to neglected dependence on the invariants, 
but rather to a weakness in the assumption that the subgrid-scale stress is solely a 
function of the velocity gradient. 

2.5 Analysis for constant coefficients 

In this section, we explore the accuracy of Eq. (11) when the coefficients are 
assumed to be constant in space. The constants are again determined through 
a least-squares procedure, this time minimizing the global error rather than the 
local error. The derivation of the global least-squares procedure is identical to that 
outlined in section 2.3, with the exception that the error is averaged over the domain 
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before it is differentiated with respect to the C,. The end result is an expression 
analogous to Eq. (12): 

Ci =< tr(m<m ; ) > -1 < tr(m,T) >, (21) 

where < > indicates a spatial average. It is important to note that when the coeffi- 
cients are determined globally, there will be non-zero error even when all five terms 
in Eq. (11) are used. 

2.5.1 Results 

As in section 2.3.1, all possible combinations of the various terms in Eq. (11) 
were tested. The correlation coefficients for the best and worst combinations of 
terms are shown in Figure 6. The actual terms corresponding to these groupings 
are listed in Table 2. 



Figure 6. Correlation coefficients for the best and worst subsets of the terms in 
Eq. (11) - constant coefficients. 

The best single term is again term 1, the Smagorinsky model. The correlation 
coefficient of the optimal groups increases as more terms are added, but the im- 
provement is rather slight (15% increase from 1 to 5 terms). The optimal correlation 
coefficients also are rather low, never exceeding 0.28. In light of the slow increase 
with additional model terms, it appears that terms 2 through 5 are not nearly as 
important as the Smagorinsky model. This conclusion can also be drawn from the 
results for the worst groupings. The correlation coefficients for the worst groupings 
are far inferior to those of the best groupings, even when 4 terms are used. The 
Smagorinsky model is the last one to be added to the worst groupings, and the 
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Number of terms 

Best combination 

Worst combination 

i 

i 

2 

2 

1,2 

2, 3 

3 

1, 2, 5 

2, 3, 4 

4 

1, 2, 4, 5 

2, 3, 4, 5 

5 

1, 2, 3, 4, 5 

1, 2, 3, 4, 5 


TABLE 2. Best and worst subsets of the model terms in Eq. (11) - constant 
coefficients. 

correlation coefficient is seen to more double when this term is added (transition 
from 4 to 5 terms). 

It is interesting to compare the results for variable and fixed coefficients (Figures 
2 and 6). Several differences are readily apparent. The level of correlation is much 
lower in the case of fixed coefficients; if only a single term is used, the correlation 
coefficient for a constant coefficient is about one half the value obtained with a 
variable coefficient. As more terms are added, the correlation improves steadily 
when the coefficients are variable, but improves little if the coefficients are con- 
stant. In the variable coefficient case, there is little to choose between the best 
and worst groupings, whereas in the constant coefficient case, the differences are 
substantial. Overall, the variable coefficient results are much better than those for 
fixed coefficients. 

The differences between the variable coefficient and constant coefficient results are 
due to differences in the scope of the minimization in the least-squares formulation. 
When the coefficients are allowed to vary in space, the error is minimized at each 
mesh point. This local minimization yields MN 3 degrees of freedom, where M is 
the number of model terms used and N 3 is the number of mesh points. When the 
coefficients are fixed in space, the error is minimized globally, and only M degrees of 
freedom are available. Evidently, the extra degrees of freedom are well utilized in the 
variable coefficient case, and superior correlations are obtained. At the same time, 
the additional degrees of freedom result in coefficients that vary greatly in space 
(recall Figure 3). This spatial dependence would be unknown in an actual LES, and 
thus the results of Figure 2 could not be realized in practice. The negative impact 
of the large coefficient variation is accounted for in Figure 6, and these results could 
be expected in practice. 

The large coefficient variation also has an important physical implication. If 
it were true that the subgrid-scale stress depended only on the velocity gradient 
tensor, then the expansion given in Eq. (11) would be complete. The coefficients 
could vary in space, but this variance would have to result from dependence on the 
invariants listed in Eq. (11). Since the coefficients are observed to vary and this 
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variation does not appear to be connected with the invariants, it must be true that 
the subgrid-scale stress at one point in space depends on more than the velocity 
gradient at the same point. While this conclusion might have been anticipated, 
the more relevant issue is to what extent the expansion in Eq. (11) captures the 
dependence of the subgrid-scale stresses on the resolved variables. In view of Figure 
6 and Table 2, it is clear that the dominant term is the Smagorinsky model. The 
remaining terms in Eq. (11) appear to be of lesser importance. In fact, if all of 
the terms are used, the correlation is only 15% higher than with the Smagorinsky 
model. Thus, at least for homogeneous isotropic flow, the expansion in Eq. (11) 
does not seem to contain much of the physical mechanisms by which the large scales 
influence the small scales. 

2.6 Summary 

A tensor relationship between subgrid-scale stress and the velocity gradient tensor 
has been developed. This relationship takes the form of a series expansion involv- 
ing products of the strain and rotation rate tensors. The expansion was used as 
a modeling hypothesis, and the latter was evaluated using direct numerical simu- 
lation data for homogeneous isotropic turbulence. The Smagorinsky model, which 
is one of the terms in the expansion, was found to be the dominant term. The 
remaining terms were found to be of lesser importance and, when included, did 
not significantly improve upon the Smagorinsky model. These results suggest that 
while the expansion is exact, the inherent assumption that the subgrid-scale stress 
depends only on the velocity gradient tensor is not well supported by the numerical 
simulation data for homogeneous isotropic turbulence at low Reynolds number . 

3. Future plans 

The conclusions drawn in the previous section apply only to homogeneous iso- 
tropic turbulence at low Reynolds number. Both the success of the model and the 
coefficient values could be Reynolds number dependent. This issue will be addressed 
by repeating the tests with higher Reynolds number DNS data. For this purpose, 
forced homogeneous isotropic simulation data is available with Reynolds number 
roughly four times greater than that used in the present study. In addition to 
Reynolds number effects, the success of the proposed model may be related to the 
flow situation. For example, it is quite possible that the model would work better in 
a shear flow where the effects of rotation are more pronounced. This possibility will 
be explored by testing the model with DNS data for homogeneous turbulent shear 
flow and for turbulent channel flow. If these results are sufficiently encouraging, 
the model will be used in an actual large eddy simulation and the results compared 
with experimental or DNS data. For the purpose of simulation, the procedure of 
§2.5 will be used to assign constant values to the expansion coefficients. 

A separate attempt will be made to use the model in conjunction with the dy- 
namic procedure of Ghosal e/ al . ( this volume). In this procedure, information con- 
tained in the resolved field will be used to estimate the value of the expansion 
coefficients as a function of space and time. This approach has the potential to re- 
cover the accuracy displayed in Figure 2 since the coefficients will be free do develop 
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any arbitrary degree of variability. 
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